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SUBSTRUCTURING  METHODS  FOR  PARABOLIC  PROBLEMS 

Maksymilian  Dryja  ' 

Abstract 

Domain  decomposition  methods  without  overlapping  for  the  approximation  of  parabolic 
problems  are  considered.  Two  kinds  of  methods  are  discussed.  In  the  first  method  systems  of 
algebraic  equations  resulting  from  the  approximation  on  each  time  level  are  solved  iteratively 
with  a  Neumann-Dirichlet  preconditioner.  The  second  method  is  direct  and  similar  to  certain 
iterative  methods  with  a  Neumann-Neumann  preconditioner.  An  analysis  of  convergence  of 
the  methods  is  presented. 
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1.  Introduction.  In  this  paper  domain  decmposition  methods  for  solving  discrete 
parabolic  problems  are  discussed.  The  discrete  problems  result  from  finite  element  and  finite 
difference  approximations  of  parabolic  problems  with  respect  to  the  space  and  time  variables, 
respectively.  For  simplicity  of  presentation  only  piecewise  linear,  continuous  approximation 
on  triangular  elements  and  the  backward  Euler  and  Crank-Nicolson  schemes  are  considered. 

Two  kinds  of  substructuring  methods,  i.e.  domain  decomposition  methods  which  do  not 
use  overlapping  subregions,  are  constructed  and  analyzed.  In  the  first  method  systems  of 
algebraic  equations  resulting  from  the  approximation  on  each  time  level  are  solved  iteratively 
by  a  conjugate  gradient  method  with  a  Neumann-Dirichlet  preconditioner.  As  in  the  elliptic 
case,  it  is  shown  that  the  rate  of  convergence  is  almost  optimal  with  respect  to  h,  the  parameter 
of  triangulation,  and  that  it  is  independent  of  the  time  step  r;  see  Section  2. 

Iterative  substructuring  methods  with  a  Neumann-Dirichlet  preconditioner  for  elliptic  finite 
element  problems  have  been  analyzed  in  [3],  [4]  and  [5]. 

The  second  method  is  direct  and  similar  to  certain  iterative  substructuring  methods  for 
elliptic  finite  element  problems  using  a  Neumann-Neumann  preconditioner;  see  [1],  [5].  A 
solution  at  each  time  level  is  obtained  in  two  fractional  steps.  The  rate  of  convergence  of  the 
method  is  of  the  order  t^'^  +  h  provided  that  r  is  proportional  to  h;  see  Sections  3  and  4. 

Most  papers  on  domain  decomposition  methods  are  devoted  to  elliptic  problems,  see  for 
example  [5]  and  the  literature  cited  therein.  Extensions  of  the  methods  to  parabolic  problems 
have  been  considered  in  [2]  and  [6]. 

2.  The  differential  and  discrete  problems.  We  consider  the  parabolic  problem.  Find 
ue  L^{0,T; H^i^))  n  C°{0,T; L'^iQ))  such  that 

(2.1a)  (^  ,<P)  +  a(u,4>)  =  ifA)  ,         (PeH^in),     a.e.      /e(o,r) 

(2.16)  {u,<t))  =  {uo,4>)  ,        (t>eL^{n) 

Here  (•,•)  is  the  Z'^-scalar  product,  a{u,v)  =  {Vu,Vv)  +  {u,v)  and  fi  is  a  bounded  polygonal 
region  in  R'^.  We  assume  that  /  G  Z^(0,T;  Z^(fi))  and  uq  £  H^{n).  The  problem  (2.1)  has  a 
unique  solution  and  is  stable,  see  for  example  [7]. 

We  solve  the  problem  (2.1)  using  a  finite  difference  method  (FDM)  and  a  finite  element 
method  (FEM)  for  the  t  and  x  variables,  respectively.  For  simplicity  only  the  backward  Euler 
and  Crank-Nicolson  and  piecewise  linear  approximations  are  discussed.  A  triangulation  of  Q 
is  constructed  as  follows:  We  first  divide  Q.  into  triangular  substructures  Q^,  with  a  parameter 
H,  and  we  then  divide  each  of  the  fi,-  into  triangles  Cj,  with  a  parameter  h.  These  form  a 
coarse  and  a  fine  triangulation  of  fi.  We  suppose  that  they  are  shape  regular  in  the  sense 
common  to  finite  element  theory.  V^{^)  is  the  space  of  continuous,  piecewise  linear  functions 
on  the  fine  triangulation,  which  are  zero  on  dO..  The  interval  [0,T]  is  partitioned  uniformly, 
t^  =  nT,n  =  0,...,N,  Nt  =  T. 

The  discrete  problem  is  of  the  form:  For  n  =  0, 1, . . . ,  A'^  —  1, 

(2.2a)  ({/^<A) +  a(f/"+',<^)  =  (/"+!, <^)  ,      <^  G  V''(fi)  , 

(2.26)  {U\(l>)  =  {uo,4>).     c^eV^^Q), 

where  [/"  =  (f/""*"^  -  U^)/t,  V^ix)  =  U{x,tn)-  The  problem  has  a  unique  solution  and  is 
stable. 


We  rewrite  (2.2)  as  the  system  of  linear  algebraic  equations:  For  n  =  0, . . . ,  TV  -  1, 

(2.3)  Air+^  =  (M  +  rA')f/"+*  =  £"+^ 

where  (Mu,  u)/?^  =  {u,v),  (A'ti,u)fi„  =  a{u,v),  u,v  £  V'^  and  m  is  the  number  of  nodal 
points  in  Cl.  It  is  easy  to  derive  a  formula  for  F""*"^  from  formula  (2.2a).  We  first  describe  an 
iterative  method  with  a  Neumann-Dirichlet  preconditioner  for  solving  (2.3). 

3.  The  Neumann-Dirichlet  iterative  method.  In  order  to  describe  the  method,  we 
first  assume  that  there  is  a  red-black  ordering  of  the  substructures  fi,.  We  will  refer  to  the  red 
and  black  subregions  as  Neumann  (N)  and  Dirichlet  (D)  substructures,  respectively.  Thus, 
no  two  substructures  of  the  same  type  share  an  edge.  For  the  case  when  there  is  no  such 
ordering,  see  the  remark  below.  Let  fiyv  and  Clp  denote  the  union  of  the  A'^-type  and  D-type 
substructures,  respectively  and  let  F  =  d^N  \  90.  and  Q.i\r  =  Q.f\f  U  F. 

We  now  stop  underlining  vectors  and  represent  the  system  (2.3)  as 


(3.1)  AU' 


^1 

0 

Al3 

0 

A22 

>i23 

Aj^ 

aL 

^33 

/?"  +  ! 


The  matrix  A\\  represents  the  couplings  between  pairs  of  degrees  of  freedom  associated  with 
the  set  fi£),  j4i3  the  couplings  between  pairs  associated  with  fi/p  and  F  and  so  on.  M  and  A' 
are  represented  in  the  same  way  i.e.  Aij  =  M,_,  +  TKij.  Let  S  be  the  Schur  complement  of  A 
with  respect  to  A33.  We  obtain 

(3.2)  S  =  5^^)  +  5(^'  , 

where  5^^)  =  A^g^  -  Al^A:^:^  A22  and  5(^'  =  /l£'  -  ^[3^1-^13.  To  define  ^^3',  we  introduce 
the  matrix  A^^^  which  corresponds  to  the  bilinear  form  {u^v)  +  Ta{u,v)  restricted  to  Qd  U  F, 
This  matrix  is  represented  as 


\    ^13       ^33     / 


where  ^l^  corresponds  to  the  couplings  between  pairs  of  variables  on  F.   /l'    '  and  /I33   are 
defined  in  a  similar  way.  We  note  that  A^^  is  the  sum  of  ^33  and  /I33  . 

Let  V^{T)  be  the  restriction  of  y''(fi)  to  F. 

Theorem  3.1.  For  all  v^  G  V^{V) 

(3.3)  (5'<^^^3,t^3)<(5t;3,^3)<7(l  +  log|^)'(5(^)r;3,t;3), 

where  7  is  a  constant  independent  of  h,  H  and  r,  (•,•)  =  {■,-)Rmi  and  mj,  is  the  number  ol 
nodal  points  on  F. 

The  left  inequality  of  (3.3)  is  obvious.  To  prove  the  right  one,  we  need  two  lemmas. 

Lemma  3.1.  For  functions  v  G  V^''(fi),  which  are  solutions  of 

(3.4)  Auvi  +  Ai3V3  =  0  ,         A22V2  +  A23V3  =  (i 
the  following  hold 

{Av,v)Rrr.   =   iSv3,V3)R'^Z    ,  iA^^HvuV3f,{Vl,V3f)fi,r,,    =   (5(^'i;3,f3)Rm3    , 


where  v  =  {vi,V2,V3)    ,  v,  £  i?"*'. 
Corollary. 

with  Vi  G  i?'"'  arbitrary. 

Lemma  3.1  and  its  Corollary  are  proved  straightforwardly  doing  some  manipulations. 
Let 

b,{u,v)  =  (w,v)L2(n,)  +  Ta^(u,v) 

where  ai(u,v)  is  the  restriction  oi  a{u,v)  to  fi,. 

Lemma  3.2.  For  the  solution  u  o/(3.4)  the  following  holds 

(3.5)  b,{u,u}  <  0(1  +  log  ^f^bj(u,u} 

where  the  summation  is  over  the  Q.j  which  are  the  N -type  neighbors  of  the  D-type  subregion 
Q,i  and  C  is  a  constant  independent  of  /i,  77,  and  t. 

Proof:  In  view  of  the  Corollary  to  Lemma  3.1,  we  have 

(3.6)  6j(u,u)  <  b,{w,w) 

where  w  is  arbitrary  at  the  nodal  points  of  $7,  and  equals  u  on  9Q,.  We  represent  w  as 
w  =  w\  ■]-  W2  -\-  w:i  where  w^  is  equal  to  w  on  side  Vij  of  fi,  including  one  of  its  end  points  x,j 
and  zero  on  50,  \  (r,j  U  {a^jj})-  Below,  we  construct  w^  such  that 

H  o 

(3.7)  b,{wj,Wj)  <  C{\  ^\og— Ybj{u,u)  ,j  =1,2,2  , 

where  bj{u,u)  is  defined  for  the  iV-type  fij  with  side  T,j,  i.e.  r,j  =  fi,  fl  fi_,.  Summing  this 
with  respect  to  j  and  using  (3.6)  we  obtain  (3.5). 

To  construct  the  w.,,  we  first  extend  the  function  u,  given  in  the  A'^-type  substructures  Qj, 
to  a  function  Wj  defined  in  a  larger  region  G  containing  the  Z?-type  substructure  fl,  and  tlie 
iV-type  substructure  flj,  such  that  w-j  =  u  in  0^,  lij  £  Hq{G)  and 

II^jIIl2(g)  <  C||w||L2(nj)  ,         \wj\h^g)  <  C'i"l//'(nj)  • 

The  function  Wj  is  not  a  finite  element  function  in  general.  Using  the  extension  theorem  from 
[8],  we  construct  a  finite  element  function  Wjh  such  that 

(3-8)  II^j/iIIl^cg)  <  C'||u||£,2(f^^)  ,     \wjk\HHG)  <  C\u\fji(Q^)  ■ 

Let  ^ij  be  a  harmonic  function  in  Q,  defined  by  its  values  on  c)Q,.  At  the  nodal  points 
of  r,j  U  {xij},  it  is  equal  to  1,  it  vanishes  at  the  remaining  nodes  and  is  defined  by  linear 
interpolation  between  the  nodes.  The  function  ^,j  €  H^^^{i},)  and  it  can  be  shown  that 

l6.lHMn.)<C(l  +  logf) 

We  now  choose  Wj  =  Ihi^tjWjh)-,  where  Ih  is  the  standard  linear  interpolation  operator. 
We  note  that,  in  view  of  (3.8), 

ll^^jllL2(n,)  <  C'll6jlli-(n.)ll^j/'llL2(n.)  <  C'll"llL2(n,)  • 


We  now  estimate  ||Vu)j||i2(nj.  Let  ^f.  be  the  average  value  of  ^,j  over  the  element  Ck-  We 
have 


and 


k  k 


in  view  of  (3.8).  Using  a  standard  error  bound  and  the  fact  that  ^,j  €  -//"'/^(fij),  we  obtain 

Summing  with  respect  to  k  and  using  (3.8)  and  a  discrete  Sobolev  inequality,  see  e.g.  [5],  we 
obtain 

E  l|V(A(6,  -  i-j)wAh(e,)  <  C(l  +  log  fnn^jhWhiu,)  <  C{1  +  log  jf\u\HHn,)  ■ 

k 


Hence 


H 


|V^.II^(Q.)^C'(l  +  log-)>|^,(„^ 


Combining  this  with  the  estimate  for  ||^j|||,2m.\)  we  obtain  (3.7). 

Proof  of  Theorem  3.1:  The  right  inequality:  It  is  enough  to  show  that 


{S^^K3,V3)  <  C{1  +  \0g^f{S^^K3,V3) 

n 


(3.9) 

Using  Lemma  3.1,  we  have 

where  vd  =  {vi,^^)'^  is  the  solution  of  (3.4)  and  vd  =  {vi,V3)'^  with  an  arbitrary  vi.  Applying 
Lemma  3.2  and  then  again  Lemma  3.1,  we  obtain  (3.9).  The  proof  is  complete. 
We  now  define  a  preconditioner  A  for  A  of  the  form 

/    ^11  0  ^13  \ 

i  =  0         ^22  A23 

\    A"^      4^       4'^^  4-  4^  A-'^  A,.    I 

\    ^13      ^23      ^33    ^  ^13"^11  '^IS    / 

Theorem  3.2.  For  any  v  £  V^{Q) 


(3.10) 


H  0  - 

{Av,v)  <  {Av,v)  <  7(1  +  log  — )  {Av,v) 


where  7  is  a  positive  constant  independent  of  k,  H  and  r. 
Proof:  It  is  easy  to  see  that 


A  =  L'DL      and      A  =  V  DL 


where 


(3.11: 


D 


(  An 

0          0      ^ 

(    ^11 

0 

0 

0 

A22       0 

, 

D  = 

0 

A22 

0 

I     0 

0     5('^'>  y 

\     » 

0 

.s- 

1  ^^ 

0 

^n  ^13  \ 

L=        0 

h 

A22  A22 

V  0 

0 

h 

The  generalized  eigenvalue  problem 

Az  =  XAz  ,  Z  =  (2l,22i  '3) 

reduces  to 

since  the  matrices  D  and  D  have  the  same  first  two  block  rows.  Applying  Theorem  3.1  we 
obtain  (3.10)  and  the  proof  is  complete. 

^From  Theorem  3.2  foUows  that  the  systems  (2.3)  with  the  matrix  A  can  be  solved  itera- 
tively  using  A  as  preconditioner  in  a  conjugate  gradient  method.  The  number  of  iterations  to 
obtain  the  solution  with  accuracy  e  is  on  the  order  of  log  -(1  +  log  y).  In  each  iteration  step 
a  system  with  the  matrix  A  is  solved.  For  that  the  factorization  A  =  {L^D)L  is  used.  In  each 
step,  this  invovles  solving  two  sets  of  the  Dirichlet  subproblems  with  the  bilinear  form  b,{u,v) 
for  the  D-type  subregions  and  a  set  of  the  Neumann  problems  for  the  A'-type  subregions. 

The  D-type  subproblems  are  independent  and  can  be  solved  in  parallel.  The  Neumann 
subproblems  are  coupled  at  the  cross  points,  i.e.  the  vertices  of  A^-type  subregions.  We  can 
use  a  block  Gaussian  elimination  or  a  capacitance  matrix  method  for  solving  the  systems 
corresponding  to  these  subproblems,  cf.  [4]. 

The  connections  at  the  cross  points  provide  a  mechanism  for  the  global  transportation  of 
information  similar  to  that  of  a  finite  element  approximation  of  an  elliptic  problem.  As  was 
shown  in  [9],  the  condition  number  of  any  method  without  such  a  mechanism  deteriorates  at 
least  as  fast  as  H~^.  For  parabolic  problems  the  situation  is  different.  As  we  will  see  below 
this  mechanism  can  be  dropped  if  r  j H^  is  bounded  cf.  [2].  In  particular,  the  set  of  Neumann 
subproblems  can  be  split  into  independent  subproblems  and  solved  in  parallel.  We  now  study 
this  case. 


Let  i4(^'  denote  the  matrix 


A22     A22, 


Note  that  for  u^v  ^  V^{^)^  which  are  the  solutions  of  (3.4), 

{S^^K^.V:,)  =  iA^^Hv2,V3f,{u2,usf  =  (1  +  t){v,u)l,(^^^  +  r(V«, V«)^2(n^,  ; 

cf.  Lemma  3.1.  To  get  the  preconditioner  A  for  A,  we  replace  /l'^^  by  A''^^  which  is  defined 
as  follows:  Let  il^^k  and  fi//  denote  the  sets  of  nodal  points  of  the  fine  mesh  in  .Qat  and  the 
coarse  mesh  in  fi,  respectively.  Let  Iqu  G  V'^{Q)  be  zero  for  x  G  ttj^i^h  \  ^H  and  equal  to  u{x) 
for  X  G  fi//.  We  introduce  a  bilinear  form  c{u,v)  by 

(3.12)  c(u,u)  =  (1 +  r)/i^     Y^    u{x)v{x)  +  T   ^   u{x)v{x) 

xenN.h  ^^"^' 

+r(V(u  -  lou)  ,  V(v  -  lov)  )^2(n^)  , 

and  the  matrix  A'^'  by 

{A^''\u2,Usf  ,(V2.Vsf)^r.,   =c{u,v), 

Let 

I(^)  =  (  ^4^     ^S  )    ,         5(^)  =  Ag^  -  Al,A^^A2s  . 

\    ^23      "^33     / 


By  using  (3.11),  we  find  that 


A  =  L'DL,         D 


(  Au       0         0 
0       A22       0 
V     0         0       5'^) 


Lemma  3.3.  For  all  v  e  V'iil),  v  =  {vi,V2,V3y 

(3.13)  foil  +  ^)~\l  +  log^)-\A^^KN,VN)  <  {A^^^vn.vn)  <  liiA^^^VN,VN) 

where  u/v  =  (^'2,^'3)'^  and  the  7,  are  constants  independent  of  h,  H  and  t. 
Proof:  We  first  prove  the  right  inequality.  It  is  easy  to  show 

Using  the  triangle  and  inverse  inequalities,  we  get 

i^v,Vv)L2^^^)  <2{\^{v  -  Iov),V{v  -  lov))^,^^^)  +  C   Yl   ^'(^)- 

xeQn 

The  right  inequality  of  (3.13)  now  follows. 
We  have 

(3.14)  (V(v  -  lov),  V(v  -  lov)  )^.(fj^)  <  2(Vt;,  Vv)^^^)  +  C   Yl   ^'(^)  • 

xe^H 


It  is  known  that 


Hl^(u.)  <  C{H-'\\v\\l^^^^  +  (1  +  log  l)\\Vv\\h(a.))  ■ 


iFrom  this  we  obtain 

(3.15)  r   Yl   ^'(^)<C'(l  +  log|^)(l  +  -^)(||t^|ll2(n^)  +  r||Vt;||i,,^^)). 

Substituting  (3.15)  into  (3.14)  and  using  the  inequality 

h'  x:   v\x)<c\\v\\l.^^^^, 

we  obtain  the  left  inequality  of  (3.13). 
Theorem  3.3.  For  all  v  G  V^((^) 

72(1  +  ;^)"'(1  +  \ogj)-\Av,v)  <  {Av,v)  <  73(1  +  \og-f{Av,v)  . 

Here  7,,  z  =  2,3  are  constants  independent  of  h,  H  and  r. 

Proof:  This  result  follows  from  Lemma  3.3  and  Theorem  3.2  and  the  representations  of 
A,  A  and  A,  see  (3.11). 

Remark.  In  the  method  discussed,  we  have  assumed  that  there  exists  a  red-black  ordering 
of  the  substructures  i.e.  no  two  substructures  of  the  same  type  share  an  edge.  To  extend  the 


method  to  the  case  when  this  condition  is  not  satisfied,  we  introduce  ND-type  substructures 
which  borders  substructures  of  both  N-  and  D-type.  In  this  case,  we  construct  a  preconditioner 
as  follows:  On  each  fi,  let  u  be  represented  as  u  =  PiU  +  HiU  where  //,u  is  the  discrete  harmonic 
extension  of  u  from  dili  in  the  sense  of  6(-,  ■)  =  (1  +  t){u,  v)  +  r(V?i,  Vv).  On  a  fi,  of  ND-lype, 
denoted  by  fi,,ArDi  we  consider  an  extension  7/,u  which  satisfies  b,^j\jD{lii,4>)  =  0  for  all  basis 
functions  <p  associated  with  the  nodal  points  x  G  ^i,ND  \  Tij  and  with  HiU  =  u  on  the  sides 
r,j  of  ^i^ND  which  are  shared  with  D-type  substructures.  On  the  remaining  sides  of  fi,,A'/:), 
we  use  homogeneous  Neumann  boundary  conditions.  Let 

d{u,v)  =  ^b,iP,u,  P,v)  +  J2f}>.N{H,u,H,v)  +  J2l}t,NDiH:U,  H,v) 

i  i  i 

and  let  5  be  a  matrix  corresponding  to  this  bilinear  form.  Here  6j,;v(-,-),  ''i,A'd(-,-)  are  the 
restrictions  of  the  bilinear  form  b{-,-)  to  substructures  of  jV-type  and  ND-type,  respective- 
ly. It  can  be  proved  that  B  is  spectrally  equivalent  to  A  except  for  a  factor  (1  -f  log^)  ; 
cf.  Theorem  3.2. 

4.  The  Neumann-Neumann  direct  method.  In  this  section,  we  discuss  a  direct 
method  for  problem  (2.1),  which  can  be  interpreted  as  a  generalization  to  the  parabolic  case  of 
a  method  with  a  Neumann-Neumann  preconditioner  previously  developed  for  elliptic  problems; 
cf.  [1],  [5].  The  problem  (2.1)  is  approximated  by  the  following  scheme. 

Forn  =  0,1,...,7V-  1, 

(4^n)  S  iUrA)^aD{{U^^"^  +  uni2,4>)^Ui^\4>),    4>eV\^) 


and 

(4.16) 


([/"+'/', <^)  +  a/v((f/"+i  +  f/"+i/2)/2,0)  =  U^^\4>) ,    <^  e  V\^) 


Here  UJ"  =  ((A"+i/2  _  U'')/t,  t/"+^/^  =  (C/"+i  -  ?7"+^/2)/r,  7"+^  =  /r+'  +/2"^^  and  /i"+'  =  0 
and  f2'^  =  0  for  X  G  ^N,h  and  x  G  ^D,h,  respectively,  a/?  and  a]si  are  the  restrictions  of  the 
bilinear  form  a{u,v)  to  Qp  and  Clj\j.  ftjsr^k  and  ^D,h  are  the  sets  of  nodal  points  in  ^i\i  and 
Q.D,  respectively,  [7"+^/^  g  V''(J7)  is  defined  by  the  values  of  [/"  at  the  nodal  points  x  G  fiyv 
and  by  the  solution  of  (4.1a)  in  CIq.  The  function  U^'^^  is  defined  similarly.  The  schemes 
(4.1a)  and  (4.1b),  taken  separately,  do  not  approximate  the  problem  (2.1).  In  spite  of  this,  we 
can  show  convergence  with  a  rate  r'/^  -|-  h  provided  that  r  is  proportional  to  h.  We  formulate 
the  stability  and  convergence  theorems  without  proofs. 

Theorem  4.1.   The  solutions  o/(4.1)  satisfy  the  following  inequality 

N-\ 
max||(y    \\i-2(Q)+  1^{\\U  +U    \\h^(Uo]  +  \W         +  <-'  \\h'(Un]) 

n=0 

<c(x:'r(ii/r^ii^+ii/r^in+iit.°in 

n=0 

where  C  is  a  constant  independent  of  h,  II  and  t. 


Theorem  4.2.   If  the  solution  u  o/(2.1)  is  sufficiently  smooth  and  t  is  proportional  to  h 
then, 

^ax  ||."||i.,„)  +  r  X:'  (||."+^/^  +  z-WUu,)  +  Ik"^^  +  ^""^'"WIhu,))  <  C{r  +  h')  , 

n  =  0 

where  z'"  =  u^  -  [/",  ^"+1/2  ^  ^n+i  _  f/n+1/2^  f;n  jg  ^/jg  solution  o/(4.1)  an  J  C  zs  a  constant 
independent  of  h,  H  and  t. 
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